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Abstract 

We  introduce  a  reduced  basis  method  that  computes  rigorous  upper  and  lower  bounds  of  the  energy 
associated  with  the  infinite-dimensional  weak  solution  of  parametrized  steady  symmetric  coercive 
partial  differential  equations  with  piecewise  polynomial  forcing  and  operators  that  admit  decom¬ 
positions  that  are  affine  in  functions  of  parameters.  The  construction  of  the  upper  bound  appeals 
to  the  standard  primal  variational  argument;  the  construction  of  the  lower  bound  appeals  to  the 
complementary  variational  principle.  We  identify  algebraic  conditions  for  the  reduced  basis  approx¬ 
imation  of  the  dual  variable  that  results  in  an  exact  satisfaction  of  the  dual  feasibility  conditions 
and  hence  a  rigorous  lower  bound.  The  formulation  permits  an  offline-online  computational  decom¬ 
position  such  that,  in  the  online  stage,  the  approximation  and  exact  certificates  can  be  evaluated  in 
complexity  independent  of  the  underlying  finite  element  discretization.  We  demonstrate  the  tech¬ 
nique  in  two  numerical  examples:  a  one- dimensional  reaction-diffusion  problem  with  a  parametrized 
diffusivity  constant;  a  planar  linear  elasticity  problem  with  a  geometry  deformation.  We  confirm 
in  both  cases  that  the  method  produces  guaranteed  upper  and  lower  bounds  of  the  energy  at  any 
parameter  value,  for  any  finite  element  discretization,  and  for  any  reduced  basis  approximation. 

Keywords:  reduced  basis,  a  posteriori  error  bound,  complementary  variational  principle,  partial 
differential  equations 


1.  Introduction 

The  theory  and  applications  of  the  certified  reduced  basis  method  —  a  model  reduction  tech¬ 
nique  that  aims  to  achieve  a  rapid  and  reliable  characterization  of  parametrized  partial  differential 
equations  —  have  advanced  considerably  in  the  past  decade  (see  Rozza  et  al.  [15],  Quarteroni  et 
al.  [13],  and  references  therein).  A  computationally  efficient  offline-online  construction  of  error 
bounds  has  been  one  of  the  main  focuses  of  the  certified  reduced  basis  method;  however,  to  our 
knowledge,  with  few  recent  exceptions  [18,  20],  the  existing  reduced  basis  error  bounds  are  with 
respect  to  some  finite  element  “truth”  which  is  assumed  to  be  sufficiently  accurate.  This  assump¬ 
tion  may  not  be  true  for  problems  with  spatial  singularities  and  in  any  event  is  often  not  verified 
in  a  rigorous  manner.  The  lack  of  reliable  feedback  on  the  validity  of  the  “truth”  can  lead  to  either 
an  inaccurate  reduced  basis  prediction  in  the  online  stage  (with  respect  to  the  infinite-dimensional 
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weak  solution)  or  overly  conservative  finite  element  “truth”  and  expensive  computation  in  the  of¬ 
fline  stage.  In  this  work,  we  introduce  a  reduced  basis  method  that  provides  rigorous  upper  and 
lower  bounds  of  the  energy  associated  with  the  infinite-dimensional  weak  solution  of  parametrized 
steady  symmetric  coercive  partial  differential  equations  under  two  assumptions:  the  equation  of 
interest  consists  of  piecewise  polynomial  forcing  functions  (or  “data”);  the  differential  operator  and 
data  admit  decompositions  that  are  affine  in  functions  of  the  parameters.  We  hence  aim  to  entirely 
remove  the  issue  of  the  “truth”  within  the  certified  reduced  basis  framework. 

Our  reduced  basis  method  appeals  to  the  complementary  variational  principle  (or  constitutive 
relation  error),  a  principle  that  has  been  successfully  used  in  the  construction  of  finite  element  error 
bounds  by,  for  instance,  Ladeveze  and  Leguillion  [7],  Ainsworth  and  Oden  [1],  and  Sauer-Budge  et 
al.  [16,  17].  More  recently,  in  the  context  of  model  reduction,  the  principle  has  been  used  in  the 
construction  of  rigorous  error  bounds  for  the  proper  generalized  decomposition  by  Ladeveze  and 
Chamoin  [6];  the  principle  has  also  been  applied  to  computational  homogenization  by  Kerfrieden  et 
al.  [5].  In  the  context  of  certified  reduced  basis  method  for  parametrized  partial  differential  equa¬ 
tions,  the  key  to  the  application  of  the  complementary  variational  principle  is  the  identification  of 
algebraic  conditions  for  the  dual  reduced  basis  space  associated  with  an  arbitrary  parameter  value. 
In  particular,  the  construction  of  the  dual  reduced  basis  space  must  be  independent  of  the  finite 
element  complexity.  We  will  demonstrate  that  this  is  indeed  possible  under  the  usual  reduced  basis 
assumption  of  the  affine  parameter  dependence. 

The  contribution  of  the  paper  is  the  reduced  basis  formulation  that  provides  upper  and  lower 
bounds  for  the  energy  associated  with  the  infinite-dimensional  weak  solution  of  steady  symmetric 
coercive  equations.  The  bounds  are  uniform,  as  opposed  to  asymptotic,  and  certifies  the  approxima¬ 
tion  for  any  parameter  value,  for  any  finite  element  resolution,  and  for  any  reduced  basis  resolution. 
In  addition,  for  our  particular  bound  construction,  we  may  associate  the  bound  gap  of  the  energy 
with  the  energy  norm  of  the  reduced  basis  solution  error,  and  hence  the  energy  bound  provides 
an  exact  certificate  of  the  solution  field.  The  method  admits  an  offline-online  computational  de¬ 
composition  such  that  the  online  computational  cost,  including  the  cost  associated  with  the  bound 
computation,  is  independent  of  the  underlying  finite  element  discretization.  As  mentioned  above, 
the  formulation  removes  the  issue  of  the  finite  element  “truth”  in  the  reduced  basis  certification 
process;  the  formulation  is  particularly  suited  for  problems  that  exhibit  spatial  singularity  in  which 
the  reliability  of  the  finite  element  “truth”  space  can  be  questionable. 

Before  we  proceed,  we  note  limitations  of  the  proposed  bound  strategy  based  on  the  comple¬ 
mentary  variational  principle;  the  first  three  are  inherited  from  the  finite  element  counterpart. 
First,  the  method  applies  only  to  symmetric  coercive  problems.  Second,  the  method  requires  the 
“data”  —  both  from  the  interior  forcing  and  boundary  conditions  —  that  is  exactly  representable 
as  piecewise  polynomials  with  respect  to  the  underlying  finite  element  triangulation.  Third,  the 
formulation  requires,  in  the  offline  stage,  a  non-standard  finite  element  approximation  of  the  dual- 
feasible  solution.  Fourth,  the  formulation  requires,  in  the  online  stage,  a  potentially  large  algebraic 
expansion,  which  could  compromise  the  online  efficiency.  We  discuss  potential  extensions  that  could 
remedy  some  of  these  limitations  at  the  conclusion  of  this  paper  in  Section  5. 

This  paper  is  organized  as  follows.  In  Section  2,  we  introduce  our  reduced  basis  method  with 
exact-solution  certificates  for  a  diffusion  equation  with  a  parametrized  diffusivity  tensor;  we  recall 
the  complementary  variational  principle,  review  a  computational  strategy  for  the  dual  variables, 
and  present  an  offline-online  computational  strategy  for  the  reduced  basis  method.  In  Section  3  we 
consider  various  extensions  of  the  method:  we  consider  a  parametrized  right-hand  side,  multiple 
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domain-dependent  parameters,  reaction-diffusion  equation,  (planar)  linear  elasticity  equations,  and 
affine  geometry  transformations.  In  Section  4,  we  demonstrate  the  method  on  two  examples:  one¬ 
dimensional  reaction-diffusion  equation  with  a  variable  diffusivity  constant;  planar-stress  linear 
elasticity  with  an  affine  geometry  transformation.  In  Section  5,  we  summarize  the  key  contributions 
and  identify  several  future  research  directions. 

2.  A  Reduced  Basis  Method  with  Exact-Solution  Certificates 

2.1.  Model  Problem:  Parametrized  Diffusion  Equation 

By  way  of  preliminaries,  we  first  introduce  a  Lipschitz  domain  14  C  with  a  boundary  par¬ 
tition  dft  =  r£)UrArforr,D  non-empty.  We  then  introduce  a  Hilbert  space  V  =  {v  £  H1(H)  : 
u|rD  =  0}  over  ft  endowed  with  an  inner  product  ( w,v)y  =  ■  Vvdx  and  the  induced  norm 

||  w || y  =  i y (w,  w)y-  We  next  introduce  a  parameter  space  V  C  Mp.  We  then  consider  the  following 
parametrized  elliptic  problem:  given  p  G  V,  find  u(p)  G  V  such  that 

a(u(p),  v]  p)  =  £{v),  \/v  G  V, 


where 


a(u>,v;p)  =  /  Vv  ■  D(p)'Vivdx,  \/w,v£V, 

Jn 

£{v)  =  /  fvdx  +  /  gvds ,  \/v  G  V. 

J  Q,  Jy  n 

for  D(p )  G  Wixd  a  symmetric  positive  definite  matrix  (that  is  invariant  over  14),  /  G  L2(ft),  and 
g  G  L2(rhv).  We  assume  that  D(p)  permits  a  decomposition  that  is  affine  in  functions  of  parameters: 
D(p)  =  Ylq= i  Qq{p)Dq  for  parameter-dependent  functions  @q:  V  — >  M  and  parameter-independent 
matrices  Dq  G  Wixd,  q  =  1, . . . ,  Q.  We  in  addition  assume  that  the  data  /  and  g  admit  piecewise 
polynomial  representations  of  a  degree  at  most  pf  over  14  and  respectively.  Recall  that  the 
solution  u(p)  G  V  is  the  infimizer  of  an  energy  functional: 

u(p)  =  arginf  Jp(w\p), 
wGV 


where 

Jp(w]  p)  =  ^ a(w ,  w ;  p)  -  £(w); 

here  the  subscript  “p”  stands  for  “primal.”  We  are  interested  in  the  energy  associated  with  the 
exact  solution  J(p)  =  Jv(u(p)',  p).  We  in  particular  consider  an  efficient  construction  of  upper  and 
lower  bounds  denoted  by  J^(p)  and  J\ y(p),  respectively,  over  the  entire  parametrer  range : 

We  will  later  see  that  the  above  bounds,  for  our  particular  construction  of  Jy!(p),  may  be  used  to 
bound  the  energy  norm  of  the  error  |||tt(/i)  —  un(p)\\\^  associated  with  a  low- dimensional  approxi¬ 
mation  un(p)  of  the  exact  solution  «(//);  here  the  energy  norm  is  defined  as  |||  •  |||M  =  p). 
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2.2.  Certified  Finite  Element  Method 
2.2.1.  Upper  Bound 

The  construction  of  an  upper  bound  of  the  exact  energy  J(p)  is  straightforward  owing  to  the 
variational  structure  of  the  elliptic  problem.  We  first  introduce  a  triangulation  Th  of  LI  and  denote 
the  skeleton  of  the  triangulation  by  dTh',  we  require  Th  to  be  chosen  such  that  f\K  £  Pp-f  (k),  Vk  £  Th, 
and  g|7  £  ¥pf  (7),  V7  £  dTh ■  Here,  PP(D)  denotes  the  space  of  polynomials  of  degree  at  most  p  over 
a  domain  D.  We  then  introduce  a  N- dimensional  subspace  Wv"  =  {u  £  V  :  v\K  £  Pp(k),  Vk  £  Th} 
on  the  triangulation  Th-  We  seek  the  finite  element  approximation:  given  p  £  T>,  find 

uM{p)  =  arginf  Jp(wAr-p)  ; 

wM<zVAr 

recall  that  the  infimizer  is  the  solution  of  the  weak  statement:  find  11:^'  (p)  £  V ^  such  that 

a(uJ^(p),  v;  p)  =  £{v)  ,  Vv  £  V ^  . 


We  now  note  that 

J+^{u)  =  =  jnf  Jp{w"-p)  >  inf  Jp(w,p)  =  J{p)  ; 

wM£yM  r  W£V 

thus,  for  any  given  p  £  V,  J+,J^(p)  is  an  upper  bound  of  the  exact  energy  J(p).  The  superscript 
“+,AA”  on  J+-^(p)  signifies  that  the  quantity  is  an  upper  bound  computed  in  a  A/’-dimensional 
finite  element  space.  We  note  that,  owing  to  the  piecewise  polynomial  assumption  on  /  and  g,  the 
interior  and  boundary  data  can  be  integrated  exactly  and  hence  we  do  not  commit  variational  crimes 
in  the  upper  bound  construction.  The  computational  complexity  of  the  upper  bound  construction 
is  0(Afk ),  where  the  exponent  k  >  1  dependents  on  the  linear  solver  strategy. 


2.2.2.  Lower  Bound:  Principle 

We  now  construct  a  lower  bound  of  the  exact  energy  J(p).  Towards  this  end,  we  first  introduce 
a  bilinear  form 

b(q,w )  =  f  q  ■  Vwdx,  Vq  £  ( L2(Ll))d ,  Vw  £  V. 

Jn 

We  then  introduce  a  broken  space  defined  on  the  triangulation  Th, 

V  =  {v  £  L2(H)  :  v\K  £  H1{k),  Vk  £  Th }  D  V.  (1) 

We  next  define  a  space  of  vector- valued  functions  Y  =  ( V)d .  We  then  introduce  a  space  of  dual 
variables  that  plays  a  key  role  in  the  certification  procedure, 

Q  =  {q  £  Y  :  b(q,w)  =  £(w),  \/w  £  H}; 

for  a  reason  that  becomes  clear  shortly,  we  refer  to  the  condition  b(q,w )  =  £(w),  Vw  £  V.  as  the 
dual  feasibility  condition.  We  finally  introduce  a  dual  energy  functional 

Jd{q;u)  =  ~7,  [  q-D~1(p)qdx,  VqeY; 

1  Jn 

here  the  subscript  “d”  stands  for  “dual.” 

We  now  state  the  key  proposition  for  our  lower  bound  construction: 
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Proposition  1.  For  any  q  £  Q  and  /jgP, 


J(p)  >  ■ 


Proof.  By  the  complementary  variational  principle, 

0<^E  J{q~  D{p)Vw)  ■  D_1(^)(g  -  D(p)Vw)dx 

q-D~1(ji)qdx+^  E  f  Vw  ■  D(p)'Vwdx  —  £/  q  ■  Vwdx 


1 

2 


KeTh 


=  -  m)  +  2a(w’ w)  -  %> w ) 


=  -  /x)  +  -a(w,  w)  -  £(w),  Vq&Q,  Vvj  £  V ; 


here  the  last  equality  follows  from  the  definition  of  the  dual  space  Q.  Since  u(p)  £  P,  it  follows 
that 

Jd{q',p)  <  ^a(u(/x),u(/x))  -  ^(«(/x))  =  J(/x),  Vg  £  Q  , 

which  is  the  desired  relationship.  □ 

The  proposition  suggests  that,  if  we  can  efficiently  construct  the  dual  space  Q,  then  we  may  use 
any  element  of  Q  to  construct  a  lower  bound.  In  addition,  as  regard  its  sharpness,  we  may  appeal 
to 

Proposition  2.  The  lower  bound  estimate  is  sharp  in  the  sense  that 

sup  Jd{q\li)  =  J(/x). 

g£Q 

Proof.  We  set  p(ji)  =  D(/i)V«(g).  We  readily  confirm  that  p(p)  £  Q  because 

b(p(p),  w)  =  b(D(p)'Vu(p),w)  =  a(u(p),w)  =  £(w),  Vw  £  V. 


In  addition, 

Jdip^fp)  =  -)-  f  Vu(p)  ■  D(p)Vu(p)dx  =  ]-  f  'Vu(p)  ■  D(n)Vu(n)dx  —  l{u(p))  =  J(p), 

and  hence  pip)  £  Q  is  the  supremizer  and  the  desired  equality  holds.  □ 

In  the  construction  of  the  reduced  basis  lower  bound  to  be  introduced  in  Section  2.3,  we  will 
only  appeal  to  the  fact  that  an  element  in  Q  that  approximates  p(p)  is  computable  in  the  offline 
stage.  In  other  words,  our  reduced  basis  lower  bound  does  not  rely  on  a  particular  finite  element 
approximation  procedure  for  p(p).  However,  in  below  for  completeness,  we  review  a  finite-element 
approximation  strategy  for  p(p)  that  exactly  satisfies  the  dual  feasibility  condition.  (See  also  Pled  et 
al.  [12]  for  other  approximation  strategies.) 
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2.2.3.  Lower  Bound:  Localization  by  Relaxation 

We  wish  to  construct  an  approximation  of  the  exact  flux  p(g)  G  Q  in  a  computationally  effi¬ 
cient  manner;  here  we  follow  the  approach  of  Sauer-Budge  et  al.  [16]  which  employs  a  two-stage 
procedure:  the  computation  of  an  approximate  inter-elemental  flux;  the  maximization  of  localized 
dual  problems  with  the  exact  dual  feasibility  constraint. 

Towards  this  end,  we  consider  an  elemental  decomposition  of  the  dual  feasibility  constraint  for 
the  space  Q: 


£(w)  —  b(q,w)=  /  fwdx+  /  gwds  —  /  q-Vwdx 

w  ri  «/  r  r-'T'  ’  ^ 

I  h 

/  (/  +  V  •  q)wdx  —  /  (nK  •  q)wds  +  (g  —  nK  ■  q)wds 

J K,  j dn\T  AT  j OkHT  at 


E 

«;eTh  L 


MweV 


Note  that  sufficient  conditions  to  satisfy  the  constraints  are,  in  the  sense  of  distribution, 

(2) 

(3) 

(4) 


— V  ■  q  =  f  in  H~\k),  V^eTh, 
hK-q  =  g  in  H~1/2(8k  nTN),  Vk  €  %, 
nK  ■  q\K  =  nKi  ■  q\K:  in  H~1/2(R  n  R1),  Vk,  k'  € 


here  R  denotes  the  closure  of  k,  and  hence  R  n  R!  is  the  face  shared  by  the  elements  k  and  n! . 
The  last  condition  requires  the  continuity  of  the  normal  fluxes  across  elemental  interfaces,  which 
introduces  a  global  coupling.  In  order  to  localize  the  problem,  we  now  introduce  a  trace  space  on 
the  skeleton  of  the  triangulation  dTi,  : 

A  =  {x  ■  x\i  e  H~1/2( 7)  ,  Vy  G  dTh\  x|7  =  9 ,  Vq  €  Tn}  ■ 


We  then  note  that,  if  we  introduce  an  arbitrary  fixed  element  %  E  A,  we  may  enforce  the  flux 
condition  on  Neumann  boundaries  (3)  and  the  normal  flux  continuity  constraint  (4)  by 

nK  •  q\K  =  crKx  in  H~l/2(dn),  Vk  G  Th  ■  (5) 


Here  <jk  is  a  function  over  each  elemental  face  of  k  G  Th  and,  for  a  face  shared  with  k'  E  7/,  and  for 
an  arbitrary  ordering  of  the  elements, 


&  K, 


—  1,  K  <  K1 , 

1 ,  otherwise 


for  a  face  on  Tjv,  oK  =  1.  In  words,  x  specifies  the  inter-elemental  flux  for  every  face  of  the 
triangulation,  which  in  turn  ensures  the  continuity  of  the  inter- elemental  normal  flux. 

We  may  express  this  form  of  the  space  Q(x)  C  Q,  characterized  by  the  interface  flux  y,  in  a 
more  convenient  form.  We  first  introduce  a  broken  bilinear  form 


q  ■  Vwdx, 


Vq  G  (L2(Ll))d,  \/vj  G  V, 
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where  we  recall  V  is  the  broken  space  defined  in  (1).  We  next  introduce  a  “jump”  bilinear  form 


aKu>xds, 


Vui  €  V,  Vx  G  A  . 


We  finally  define 

Q(x)  =  {q£Y  :  b{q ,  w)  =  £(w)  +  c(w,  x),  Vw  G  V}  . 


Note  that,  in  order  to  ensure  that  the  space  is  not  empty,  the  interface  flux  x  £  A  must  be 
equilibrating  in  the  sense  that  0  =  £(1K)  +  c(lK,x),  Vk  G  Th,  where  1K  =  1  on  k  and  1K  =  0  on 
fi\  k.  We  emphasize  that  Q(x)  C  Q  for  any  equilibrating  flux  x  £  A. 


2.2.4 ■  Lower  Bound:  Computation 

We  now  consider  a  finite  element  approximation  of  the  interface  flux  X(p)  G  A  and  the  dual 
variable  p(p)  G  Q(X(p)).  Towards  this  end,  we  first  introduce  a  finite  element  flux  space  =  {A  G 
A  :  A|7  G  Pp(7),  V7  G  dTh}-  (Note  that  the  dimension  of  this  space  is  not  J\f  but  is  of  0(N).)  We 
then  seek  an  approximate  interface  flux  associated  with  the  finite  element  solution  \p)  G  : 
hnd  X^ (p)  G  A  such  that 

c(w,  X^  (p))  =  a(w^ (n),w]  p)  —  £{w)  ,  Vw  G  V ^  .  (6) 

We  solve  this  so  called  equilibration  problem  using  the  method  of  Ladeveze  and  Leguillon  [7]  and 
in  particular  its  high-order  extension  by  Ainsworth  and  Oden  [1],  which  has  a  computational  cost 
that  scales  linearly  with  the  number  of  finite  element  unknowns. 

We  next  introduce  a  broken  finite  element  space  of  interior  fluxes:  =  ( 'yY)d .  We  now  wish 

to  construct  a  dual-feasible  subspace  (\^ (p))  C  Q(X^ (p))  C  Q, 

QAf(XAr(p))  =  {q  G  YN  :  b(q,  w )  =  £(w)  +  c(w,  X N  (p)),  \/u>  G  V}  .  (7) 

We  must  satisfy  exactly  the  dual  feasibility  condition  (2),  (3),  and  (4).  We  also  recall  that  we  may 
replace  the  constraints  (3)  and  (4)  by  a  single  condition  (5).  Hence,  the  sufficient  condition  for 
dual  feasibility  in  strong  form  is 

— V  ■  q  =  f  in  #-1(ac),  Vk  G  Tft  , 

nK-  q  =  c TKXjV(p)  in  H~1/2(dn),  Vk  G  %  ■ 

Note  that  our  volume  data  f\K  G  PPp(k)  and  boundary  data  XJ^(p)\^  G  Pp(7),  7  G  dn.  In  addition, 
by  our  choice  of  the  space  for  q,  V  •  q\K  G  Pp_1(k).  It  follows  that  we  need  only  test  the  equation 
over  k  against  the  finite-dimensional  space  Pp_1(k)  (forp  >  pf  + 1)  and  not  the  infinite- dimensional 
space  H1{k).  Similarly,  we  need  only  test  the  equations  over  dn  against  finite- dimensional  spaces 
IP75 (7) ,  7  G  dn.  Thus,  we  can  construct  —  implicitly  using  a  finite  number  of  constraints  —  the 
space  QJ^(XY'[p))  c  Q  that  satisfies  exactly  the  dual  feasibility  condition.  The  existence  of  at 
least  one  q  G  Y^  that  satisfies  the  dual  feasibility  conditions,  and  hence  the  non-emptiness  of 
Q^(X ^{p)),  is  discussed  in  [16]. 

We  may  now  readily  construct  a  finite  element  approximation  of  a  lower  bound  of  the  energy. 
We  first  solve  (6)  for  the  approximate  flux  function  A ^ (p)  associated  with  the  finite  element  solution 
u-^(p).  We  then  compute  the  finite  element  approximation  of  the  dual  variable 

A)  =  arg  sup 

q"  S  Q"  (AAr  (m)) 
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We  finally  evaluate  the  associated  lower  bound 

=  Jd{pM(p)-,p)  <  J(p); 

the  inequality  is  a  direct  consequence  of  Proposition  1  and  {p))  C  Q.  The  superscript 

,J\f”  on  J~,j^(p)  signifies  that  the  quantity  is  a  lower  bound  computed  in  a  0 (AO-dimensional 
finite  element  space.  Note  that  the  dual  maximization  problem  is  a  quadratic  program  with  linear 
constraints,  which  reduces  to  element-wise  saddle  problems  that  can  be  solved  efficiently.  The 
computational  complexities  of  the  inter-elemental  flux  calculation  and  local  dual  problems  are  of 
0{N)  (provided  that  the  primal  solution  it^(p),  which  is  used  in  the  calculation  of  X^(p),  has 
already  been  computed). 

2.3.  Reduced  Basis  Method 
2.3.1.  Upper  Bound 

We  now  introduce  a  reduced  basis  method  that  provides  rigorous  upper  and  lower  bounds 
of  the  exact  energy.  As  before,  the  construction  of  an  upper  bound  is  straightforward  owing  to 
the  variational  structure  of  the  elliptic  equation.  We  first  introduce  a  hierarchical  A-dimensional 
reduced  basis  space 

Vn  =  span{uA/'(^(n)),  p(n)  G  MN}  =  span{^}^=1  C  VM , 

where  Mn  =  {p(n)}n=i  is  the  set  of  N  reduced  basis  parameter  points  (for  formally  N  <  Af 
but  typically  N  <C  A f),  and  is  the  orthonormal  basis  with  respect  to,  say,  Our 

reduced-basis  approximation  is  given  by 

uN{p)  =  arginf  Jp{wN\  p)  ; 

wn£Vn 

again,  the  infimizer  is  given  by  the  Galerkin  statement:  find  un{p)  G  Vn  such  that 

a(uN(p),vN;p)  =  £{vn)  ,  Vu/v  £  Vn  ■ 

Our  reduced  basis  upper  bound  is  then  given  by 

JM  =  Jp(uN(p)-,P )  =  inf  Jp(wN;p )  >  inf  Jp(ioM-p)  =  J+M {p)  >  J(p)  ; 

wNevN  wN  evN 

The  superscript  and  subscript  on  J^(p)  signifies  that  the  quantity  is  an  upper  bound  computed 
in  a  Ar-dimensional  reduced  basis  space.  The  reduced  basis  upper  bound  is  looser  than  the  finite 
element  upper  bound. 

The  computation  of  the  reduced  basis  upper  bounds  permits  the  standard  offline-online  com¬ 
putational  decomposition.  (See,  for  example,  Rozza  et  al.  [15].)  In  the  offline  stage,  we  first  solve 
N  finite  element  problems  to  obtain  basis  functions  u-^(p^),  n  =  1, . . . ,  N,  for  the  reduced  basis 
space.  We  then  orthonormalize  the  functions  to  obtain  ,  n  =  1  ,...,N.  We  next  appeal  to 
the  decomposition  of  the  diffusion  coefficient  that  is  affine  in  functions  of  parameter  and  compute 
parameter-independent  matrices  Ay  G  M.NxN ,  i,j  =  1 ,...  ,d,  and  vector  F  €  with  entries 

(A ij)mn=  [  m,n  =  l,...,N  , 

Jn  oxi  oxj 

F m  =  t(€m),  m  =  l,...,N  . 
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In  the  online  stage,  we  take  a  three-step  procedure.  We  first  compute  a  parameter-dependent 
matrix 


d 

Mp)  =  /;o'{/'>Ad- 

ij= 1 

We  then  solve  a  IV  x  IV  reduced  basis  problem:  find  oi*N{p)  £  such  that 

A  (n)a*N(n)  =  F 

We  hnally  evaluate  the  reduced  basis  upper  bound 

Jn(v)  =  ~\(aN(p))T  A(p)a*N(p)  • 

The  computational  complexity  of  the  online  stage  is  0(d‘2N2)  +  0(N 3)  and  is  in  particular  in¬ 
dependent  of  the  finite  element  complexity  of  0(J\f)  (or  the  exact  infinite-dimensional  problem 
complexity  of  infinity).1 

2.3.2.  Lower  Bound 

We  now  construct  a  reduced  basis  approximation  of  the  energy  lower  bound.  We  first  define 
the  reduced  basis  space  associated  with  the  dual  variable 

Yn  =  span {pA/'(M(n))j  ^(n)  €  mn}  =  span {r^}^=l  C  YM . 

Here,  { r }%=1  is  the  orthonormal  basis  with  respect  to,  say,  (•,  -)(L2(n))d;  we  in  addition  introduce 
the  change  of  basis  matrix  tc  €  WNxN  that  satisfies 

N 

rtf  =  X]  ^  {P(n'))^n'n,  n=l,...,N. 

n'= 1 

We  then  define  the  key  space  for  the  lower  bound  computation:  the  reduced  basis  dual  space  with 
the  dual  feasibility  constraint, 

Qn  =  {q  G  Yn  :  b(q,w)  =  £(w),  \/w  £  V}  ;  (8) 

we  again  must  satisfy  exactly  the  dual  feasibility  conditions.  Towards  this  end,  we  substitute  the  re¬ 
duced  basis  approximation  of  the  dual  variable  pjy  =  J2n-i  rtf  PNn  =  Yln= l  Yln’= l  (PM^n’nPNn 
to  express  the  dual  feasibility  condition  (8)  in  terms  of  the  reduced  basis  coefficients  (3n  £  M N : 

N  N 

b(pN,w)  =  b(^2  ^2  PM (Piri^Un'nPNn,  w)  =  i{w),  Vw  £  V.  (9) 

n=  1  n'= 1 


1The  0(d2N2)  term  associated  with  the  assembly  process  is  an  upper  bound  of  the  0(QaN2)  term  associated  with 
the  “standard”  expression  for  the  reduced-basis  assembly,  where  Qa  is  the  number  of  terms  in  the  affine  expansion 
of  a(-,  •;  •)  [15],  for  the  spatially  invariant  diffusivity  tensor  D(y)  £  Rdx  . 
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On  the  other  hand,  we  appeal  to  the  bilinearity  of  the  form  &(•,  •)  and  the  fact  that  each  p^(n(n'))  G 
Q  —  which  implies  b(p^ \pLtn')),w)  =  t(w),  Vw  G  V,  n!  =  1, . . . ,  N  —  to  obtain 

N  N  N  N 

b(£2  E  PM^(n'))^n'nPNn,w)  =  E  E  Un'nPNnHP^ {H{n')) ,  w) 

n=l  n'=l  m=l  n'=l 

N  N 

-££<*  nPNnZ{w),  VtC  G  V.  (10) 

77.— 1  n'=l 

The  comparison  of  (9)  and  (10)  shows  that  the  sufficient  condition  for  p n  G  Qn  is 

N  N 

^  ^  ^  ^  W n'n^Nn  =  1  • 
n=l  n'=l 

With  the  constraint,  the  effective  dimension  of  the  reduced  basis  coefficient  space  is  N  —  1.  We 
now  define  our  reduced  basis  energy  lower  bound  as 

j^(p)=  sup  Jd{qN\p)=  sup  Jd(qN ;m); 

<JjvSQjv  qj veYjv 

£n  =  l  £n'  =  l  UJn'n0Nn  =  1 

The  superscript  and  subscript  on  J^(/r)  signifies  that  the  quantity  is  a  lower  bound  computed 
in  a  IV-dimensional  reduced  basis  space.  We  note  that  the  reduced  basis  lower  bound,  unlike  the 
upper  bound,  can  be  tighter  than  the  finite  element  lower  bound.  This  is  because  the  function 
Pn{p)  is  optimized  for  the  tightest  bound  from  a  set  of  dual  feasible  functions  in  Qn  C  Q  whereas 
the  finite  element  dual  variable  p^(n)  —  computed  by  the  procedure  described  in  Section  2.2.4  — 
is  chosen  from  a  subspace  ((w  (p))  C  Q  restricted  by  the  choice  of  W  (/r). 

The  evaluation  of  the  reduced  basis  lower  bound  also  permits  an  offline-online  computational 
decomposition.  In  the  offline  stage,  we  first  solve  N  finite  element  dual  problems  to  obtain  basis 
functions  jA(/t(n/)),  n'  =  1, . . . ,  N,  for  the  reduced  basis  space  Tv-  We  then  orthonormalize  the 
functions  to  obtain  rrff ,  n  =  1, . . . ,  IV,  and  the  change  of  basis  matrix  co  G  M.NxN .  We  next  sum 
the  leading  index  of  the  change  of  basis  matrix  to  form  G  with  entries 

N 

wnUm  =  E  Wn'n’  U  =  -1’  ’  '  •  ’  N- 

n'=l 

We  then  compute  parameter-independent  matrices  K,j  G  M>NxN ,  i ,  j  =  1 ,d,  with  entries 

(K ij)mn  =  j  Vm,irfn,jdx  , 

where  77E  denotes  the  i-th  vector  component  of  the  m-th  basis  function. 

In  the  online  stage,  we  take  a  three-step  procedure.  We  first  form  a  parameter-dependent  matrix 
K(/x)  G  RNxN  defined  by 

d 

k(a0  =  E 

i,j= i 
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We  then  solve  a  N  x  N  quadratic  program  with  a  single  linear  constraint 

Pn(j*)  =  argsup  K(ji)PN  , 

pNmN 

(wsum)TAv= 1 

which  requires  the  solution  of  a  ( N  +  1)  x  ( N  +  1)  saddle  system.  We  finally  evaluate  the  lower 
bound 

The  computational  complexity  of  the  online  stage  is  0{d2N2)  +  0((N  +  l)3)  and  is  in  particular 
independent  of  the  finite  element  complexity  of  0(J\f)  (or  the  exact  infinite-dimensional  complexity 
of  infinity). 

Remark  1.  Here  we  choose  the  same  set  of  parameter  values  for  the  primal  snapshots  that  spans 
Vjv  and  the  dual  snapshots  that  spans  Y,v-  In  general  the  parameter  values  associated  with  the 
primal  and  dual  snapshots  need  not  be  the  same.  In  particular,  for  online  efficiency  —  that  is 
to  reduce  N  —  it  is  in  general  advantageous  to  consider  different  parameter  values  because  the 
parametric  manifold  associated  with  the  primal  and  dual  solutions  can  be  quite  different.  However, 
for  offline  efficiency,  it  is  advantageous  to  consider  the  same  snapshot  parameter  values  because 
the  computation  of  the  finite  element  dual  solution  by  the  equilibration  procedure  requires  the 
primal  solution.  We  recall  that  the  complexity  of  the  primal  solve  is  0(Afk),  k  >  1,  where  as  the 
complexity  of  the  equilibration  procedure  is  0(J\f). 

2. 3. 3.  Energy  Norm  of  the  Error 

The  upper  and  lower  bounds  of  the  energy  functional  may  be  used  for  the  global  certification 
of  the  reduced-basis  solution  in  the  energy  norm  |||ie|||M  =  y7 a(w,  w,  //): 

Proposition  3.  The  upper  and  lower  bound  of  the  energy  functional  is  related  to  the  energy  norm 
of  the  error  in  the  sense  that,  for  any  [i  6  V, 

1 1  Km)  —  'uw(m)III  h  ~  ~  •Ijv(m))  • 

Proof.  We  appeal  to  the  definition  of  the  energy  norm,  Galerkin  orthogonality,  and  the  definition 
of  the  energy  functional  to  obtain 

1 1  Km)  -  un(h)  lll^  =  a(u(n)  -  uN(n),u{n)  -  uN(n)-,n)  =  a{u{n),u{nY,n)  -  a(uN(n),  n) 

=  —2 fj,)  +  2 /x)  <  2 ( —  J N(n))  ; 

this  concludes  the  proof.  □ 

3.  Extensions 

3.1.  Extension  1:  Parametrized  Right-Hand  Side 

We  now  consider  an  extension  of  the  method  for  the  case  in  which  the  right-hand  side  data  £(■) 
is  parametrized  but  permits  a  decomposition  that  is  affine  in  functions  of  the  parameter: 

s 

m)  =  ®s{r)£s(v)  (11) 

S=1 
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for  parameter-dependent  functions  ©s  :  V  — >  M  and  parameter-independent  functionals  £s  G  V1 , 
s  =  1, . . . ,  S.  The  primal  formulation  requires  no  modifications,  and  the  construction  of  the  upper 
bound  follows  from  the  same  variational  argument  as  before. 

To  construct  a  lower  bound,  we  appeal  to  the  linearity  of  the  dual  solution  p(p)  on  the  right 
hand  side.  We  first  note  that  the  dual  feasibility  space  is  now  parameter-dependent  and  is  given  by 

s 

Qn(p)  =  {q£YN  :  b(q ,  w)  =  ^  @s(p)£s(w),  \/w  G  V}.  (12) 

S=1 

We  next  substitute  the  form  of  our  reduced  basis  approximation 

Pn  =  ' @Nn  =  J2n=i  (.V(n'))un'n/3Nn  to  the  expression  for  the  constraint,  ap¬ 

peal  to  the  bilinearity  of  b(-,  •),  and  invoke  p^  (p(n'))  £  Q{p{n'))i  n  =  1 ,...  ,N  —  which  implies 
bipAr(p-(nl)),w)  =  Yls= i  ©s(m)4(w),  Viu  G  V,  n'  =  1, . . .,  N  —  to  obtain 
N  N  N  N 

«EE*"  EE“  n'nPNnKp^ \p(n'))  ,  w) 

n= 1  n'=l  n=l  n'=l 

S  N  N 

=  EEE  ©s  iP(n')  ^^n'n/^Nn^-s  (^)  5  ViU  G  Id 

s=l  n=  1  n'=l 

The  comparison  of  (12)  and  (14)  shows  that  a  sufficient  condition  for  p n  G  Qn(p)  is 

N  N 

^  ^  ^  ^  ®s(p(n'))^n'n/3Nn  =  ©s(m) >  S  =  1, .  .  .  ,  iS  ; 

n=l  ?i'=l 

the  set  of  5  constraints  is  a  generalization  of  the  single  constraint  Yln=i  Yln’=i  un'nPNn  =  1  for 
the  non-parametrized  right-hand  side  case.  Note  that  the  effective  dimension  of  the  reduced  basis 
space  is  N  —  S  (assuming  the  S  constraints  are  linearly  independent);  we  consequently  require  that 
N  >  S. 

The  offline-online  computational  decomposition  follows  from  the  non-parametrized  right-hand 
side  case  with  one  exception:  in  the  online  stage,  we  simply  impose  S  constraints  instead  of  the 
single  constraint.  The  computational  complexity  of  the  online  stage  is  0(d2N2)  +  0((N  +  S')3), 
where  N  +  S  is  the  size  of  the  saddle  system  with  the  S  constraints.  Note  that,  as  the  effective 
dimension  of  the  reduced  basis  space  is  N  —  S,  the  dimension  N  may  need  to  be  larger  than  the 
non-parametrized  right-hand  side  case. 

Remark  2.  Instead  of  the  procedure  described  above,  we  may  construct  S  parameter- independent 
dual  feasibility  spaces,  each  corresponding  to  the  specific  £s.  Specifically,  we  may  construct  separate 
spaces  QSjn  =  {q  G  YSi N  :  b(q ,  w)  =  ts (w),  \/w  G  V},  s  =  1, . . . ,  S,  with  YS)N  =  span{#f  (M(n))ln=i 
associated  with  is. 


(13) 

(14) 


3.2.  Extension  2:  Multiple  Domains 

We  now  consider  an  extension  of  the  method  for  the  case  in  which  the  bilinear  form  o(-,-) 
constitutes  of  contributions  from  Jidom  subdomains,  D,(k\  k  =  1 , ... . .  /Tdond 


a(w,  v ; 


Vu  •  D{k\p)Vwdx, 


Vw,  vGf; 
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we  assume  symmetric  positive  definite  matrices  D^k\n)  £  M.dxd,  k  =  1, . . . ,  i^dom;  are  constant  over 
each  subdomain  (but  in  general  different  across  subdomains).  The  construction  of  the  upper  bound 
follows  from  the  same  variational  argument,  and  the  primal  formulation  requires  no  modifications. 
To  construct  a  lower  bound,  the  dual  energy  functional  is  modified  to 

Jd(q;n)  =  ~lY^  [  q- (D{k)(fj)ylqdx  . 

2  J^k) 

With  the  redefinition,  we  may  follow  the  proof  of  Proposition  1  and  still  show  that 

Jd{q-,v)  <  Vq^Q, 


where  Q  =  {q  £  Y  :  b(q,w)  =  £(w),  \/w  £  V}  as  before.  We  in  particular  note  that  Q  is 
independent  of  the  diffusion  coefficients  k  =  1, . . .  ,K^om.  It  follows  that,  for  any  given 

^  £  V,  we  may  find  a  dual  feasible  function  £  Q  using  the  finite  element  procedure  described 

in  Section  2.2.4.  In  addition,  the  constraints  for  the  reduced  basis  coefficients  are  unchanged: 

En=l  E^=l  Un'nPNn  =  1- 

We  need  to  make  only  minor  modifications  to  the  offline-online  computational  procedure.  In 
the  offline  stage,  we  first  solve  N  finite  element  dual  problems  to  obtain  snapshots  p^(f U(n)),  n  = 
1, . . . ,  N.  We  then  orthonormalize  the  functions  to  obtain  the  basis  {rrff  }‘Ei  for  the  reduced  basis 
space  Yn,  and  the  change  of  basis  matrix  u  £  WNxN .  We  next  compute  parameter-independent 
matrices  K?E  £  lS&NxN ,  i,j  =  l,...,d,k=l,...,  Kdom ,  with  entries 


(K 


Wn  _ 


V 


law 


r,'T  .  /7'T  .  //  a’ 


In  the  online  stage,  we  first  form  a  parameter-dependent  matrix  K(/i)  £  M.NxN  defined  by 


KM=  £  ^(cW(c)).yK<‘>. 

fc=l  i,j=  1 

We  then  solve  the  N  x  N  quadratic  program  with  the  single  linear  constraint  and  form  the  lower 
bound  as  before.  The  computational  complexity  in  the  online  stage  is  0(Kdomd2N2)  +  0((N  +  1)3). 


3.3.  Extension  3:  Reaction- Diffusion  Equation 

We  now  consider  an  extension  of  the  method  to  the  reaction-diffusion  equation.  The  bilinear 
form  associated  with  the  equation  is 

a(w,v;ix)=  /  (Vv  ■  D(n)'Vw  +  c(/-i)vw)  dx 
Ja 

and  the  associated  energy  functional  is  Jv{w\  //)  =  ^a(w,w,  /a)  —  £(w).  The  solution  u(//)  £  V 
such  that  a(u,v;/i )  =  £(v),  Vu  £  V,  is  the  inhmizer  of  the  energy  functional.  Thus,  owing  to  the 
variational  structure,  the  construction  of  an  upper  bound  requires  no  modifications. 

To  construct  a  lower  bound,  we  first  introduce  the  second  dual  variable  z  £  V,  where  we  recall 
V  is  the  broken  space  defined  in  (1).  We  then  redefine  the  bilinear  form  that  induces  the  dual 
feasibility  constraint: 

b((q,  z),w )  =  /  (q  ■  Vw  +  zvj)  dx,  Vq  £  Y,  \/z  £  V,  \/w  £  V  . 

Ja 
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We  next  redefine  the  dual  feasible  space 

Q  =  {(q,  z)  £  Y  xV  :  b((q,  z ),  w)  =  £{w),  Vru  6  b}  . 

We  finally  introduce  the  dual  energy  functional 

Jd((q,  z);  fj>)  =  -X  [  (<?  ■  D~l(p)q  +  c^^zz)  dx,  Vq  6  7,  Vz  eV  . 

1  Jn 

We  may  then  show 

</(X)  >  Jd((g,  z);  n),  \/(q,  z)  eQ  , 

following  the  same  argument  as  the  proof  of  Proposition  1. 

In  order  to  construct  a  finite  element  approximation  of  the  dual  variables  (q,  z)  £  Q,  we  follow 
the  same  localization  procedure.  (See  Sauer-Budge  and  Peraire  [17]  for  details.)  Namely,  we 
first  compute  the  equilibrating  inter  elemental  flux  (p)  £  A  as  defined  by  (7)  (for  the  a(-,  •;  p) 
associated  with  the  reaction-diffusion  equation).  We  then  identify  the  strong  form  of  the  local 
constraints 

—V-q  +  z  =  f  in  77_1(k),  Vk  £  %  , 

nK  ■  q  =  crK\M{n)  in  H~1/2(dn),  \/k  £  %  ■ 

Here,  for  f\K  £  P P/(k),  A-,v”(/u)|7  £  Pp(y),  Vy  £  dn,  and  p  >  pf ,  we  may  choose  q\K  £  Pp(k)  and 
z\K  £  Pp_1(/{);  we  then  realize,  as  before,  we  need  only  test  the  equation  over  k  against  Pp_1(k)  and 
the  equations  over  dn  against  Pp(y),  7  £  Ok.  Thus,  we  can  express  the  dual  feasibility  condition  as  a 
finite  dimensional  constraints.  The  existence  of  at  least  one  solution  ( q ,  z)  £  Yv  x  that  satisfies 
the  dual  feasibility  conditions,  and  hence  the  non-emptiness  of  QJ^(XJ^(p)),  is  a  consequence  of  the 
existence  condition  for  the  diffusion  equation;  the  presence  of  the  dual  variable  associated  with 
the  reaction  term,  z  £  ,  increases  the  number  of  unknowns  while  the  number  of  constrains  is 

unchanged.  We  then  seek  {p^{p),  rM  {p))  =  argsup, M (P))  <?d((q^ ,  z^)\  p). 

We  now  deduce  the  constraints  in  the  reduced  basis  setting.  We  first  introduce  an  orthonor- 
malized  basis  of  the  reduced  basis  space  Y/v,  {(77^,  (£{  )}^Li!  the  change  of  basis  matrix  u  £  RNxN 
satisfies 

/  N  N 

(V^iCn)  =  I  ^2  Y  ^ ^{n'))un'n 

Vn'=l  n'= 1 

We  express  our  reduced  basis  approximation  as  pn  =  J2n= 1  rfn  Pnu  =  J2n= 1  1 {^{n'^n'n^Nn 

and  rN  =  J2n=iC n  Pnu  =  En=i  Yjn'=i  {V{n'))un'nPNn,  appeal  to  the  bilinearity  of  &(•,■), 
and  (p* {p(n')),'r*r (V(n')))  £  Q  —  which  implies  b((p^  {p^),  r*  (p^)) ,w)  =  £(w),  \/w  £  V, 
n!  =  1, . . . ,  N  —  to  obtain 

N  N  N  N 

b{(pN,rN),w )  =  b((Y  Y  EE^  (/^■(n')  'j^ri'n^Nn)  j  ^0 

n=ln'=l  n=  ln'=l 

N  N 

=  EE  ^n'npNnb{(p^{p(n')),1M{p(n'))),w) 

n=l  n’= 1 
N  N 

—  ^  ^  "  Mn'nfiNn^ipu)  i  £  V  , 

»l=l  7l'  =  l 
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we  recognize  the  sufficient  condition  for  pn(p)  £  Q  is  JJn=iUn'=i  ^n'nPNn  =  1  —  the  same 
condition  as  the  diffusion-only  case. 

The  offline-online  computational  decomposition  is  similar  to  the  diffusion-only  case.  In  the 
offline  stage,  we  first  solve  N  finite  element  dual  problems  to  obtain  basis  functions  p^{p(n')) 
and  n'  =  1  We  then  orthonormalize  the  functions  to  obtain  vtf  and  ,  n  = 

1, . . . ,  N,  and  the  change  of  basis  matrix  oj  G  M.NxN .  We  next  sum  the  leading  index  of  the  change 
of  basis  matrix  to  form  w®um  =  Y2n'=i  urin-  We  then  compute  parameter- independent  matrices 
K ij  G  WNxN,  iff  =  1, . . . ,  d  and  C  G  M,  with  entries 

(Kjj)rnn  =  /  'llm,irtn,jdx  and  Cmn  =  /  CmCrf dx  . 

J  J  f2 

In  the  online  stage,  we  take  a  three-step  procedure  as  before.  We  first  form  parameter-dependent 
matrices  K(/j)  G  M.NxN  and  C(p)  G  WNxN  defined  by 

d 

K(h)  =  and  C(/i)  =  c_1(/i) C  . 

i,j=  1 

We  then  solve  a  N  x  N  quadratic  program  with  a  single  linear  constraint 

Pn(v)=  argsup  +  C(p))/3N  , 

/3N£RN  1 
(ojsum)T/3iv  =  l 

which  requires  the  solution  of  a  [N  +  1)  x  (N  +  1)  saddle  system.  We  finally  evaluate  the  lower 
bound 

W  =  +  C Gu))/4(m)  • 

The  computational  complexity  in  the  online  stage  is  0(d2N 2)  +  0({N  +  l)3). 

Remark  3.  We  have  here  used  the  same  reduced  basis  coefficients,  /3jv,  for  pjy  =  rfn  Pnu  and 

rN  =  En= 1  /3jvn-  We  may  instead  consider  different  coefficients  for  the  two  spaces.  The  coupled 
formulation,  as  considered  in  this  work,  is  typically  more  online  efficient  as  p(p)  and  r(p)  tend 
to  be  correlated.  The  decoupled  formulation  is  less  online  efficient  but  is  more  offline  efficient:  it 
produces  an  online  system  of  size  2N,  but  the  reduced  basis  space  has  a  larger  dimension  for  a 
given  number  of  offline  snapshots. 

3-4 ■  Extension  4-  Affine  Geometry  Transformation  (Reaction- Diffusion) 

We  now  consider  an  application  of  the  method  to  problems  with  an  affine  geometry  transfor¬ 
mation.  By  way  of  preliminaries,  we  introduce  a  parameter-independent  reference  domain  If  and  a 
parameter-dependent  transformed  domain  D(p).  A  point  x  G  Cl  is  mapped  to  a  point  x  G  D(p)  by 
an  affine  transformation 


x  =  G(x\  ff)  =  T(p)x  +  x0(p), 

where  T(p)  G  Mdxd  is  the  Jacobian  of  the  transformation,  and  xo(p)  G  is  associated  with  the 
translation. 
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We  may  readily  transform  the  forms  associated  with  the  parameter-dependent  domain  Q(/j)  to 
the  forms  associated  with  the  parameter-independent  reference  domain  Q.  (We  refer  to  Rozza  et 
al.  [15]  for  more  detailed  discussion  in  the  standard  reduced  basis  context.)  We  first  transform  the 
bilinear  form  a(-,  ■;  g):  for  w  =  w  o  G(-;  g)  and  v  =  v  o  £?(•;  g), 


a(w,  v: 


IQ  L 


[Vu  •  D(g)Vw  +  c(g)vw]  dx  =  /  I'Vv  ■  D(fi)'Vw  +  c(fj,)wv  dx  =  a(w,v\ g), 


where  D(g)  =  det (T(g))T  T (g)D(g)T  1(g)  and  c(g)  =  det {T(g))c(g).  We  then  transform  the 
linear  form  £(■;  g):  for  v  =  v  o  G(-;  g), 


£{v\g)  =  /  f(g)vdx  +  /  g(g)vds  =  f(g)vdx  +  g(g)vds  =  I(v;  g), 

J  r2(^.)  */ r  jy  (/z)  J  f  jy 

where  f(g)  =  (det (T(g))f(g))  o  G{-;g)  and  g(/u)  =  (det(TrN(g))g(g))  °  G{-;g).  We  in  addition 
define  the  primal  energy  functional  on  the  reference  domain:  Jp{-;g)  =  \a{-,-;g)  —  £{-;g).  The 
transformation  for  o(-,  •;  g)  and  £(•;  g)  are  in  fact  the  standard  transformation  in  the  reduced  basis 
context  [15]. 

We  now  introduce  a  key  transformation  for  the  lower  bound  construction.  We  transform  the 
bilinear  form  b(-,-;g),  which  is  now  parameter  dependent  due  to  the  presence  of  0(/x) ,  as  follows: 
for  q  =  (det (T(g))T~T(g)q)  o  G(-;  g),  z  =  (det (T(g))z)  o  G(-;  g),  and  w  =  w  o  G( ■;  g), 

b((q,  z),w;  g)  =  /  (q  ■  Vw  +  zw)  dx  =  /  (q  ■  T~1(g)Vwdet(T(g))  +  zwdet(T(g))\  dx 
JQ{n)  J Q  '  ' 

q  ■  Vw  +  zw'j  dx  =  b((q,  z),w). 

We  note  that  the  transformed  bilinear  form  b(-,  •)  is  in  fact  independent  of  the  parameter  g.  Ac¬ 
cordingly,  we  may  identify  the  space  of  dual  feasible  functions  in  the  reference  domain 

Q(Cl;g)  =  {( q,z )  G  Y(Ct)  x  V(Cl)  :  b((q,z),w)  =  £(w;g),  \/w  G  R(f7)},  (15) 

where  the  parameter-dependence  of  the  dual-feasibility  condition  arises  only  through  the  right- 
hand  side,  £,  and  can  be  treated  using  the  technique  in  Section  3.1.  We  readily  confirm  that,  if 
{q,z)  G  Q{Cl;g),  then  (q,z)  =  (((det(T(^)))"iTT(/r)g)  o  G_1(-,  g),  ((det(T(/i)))_15)  o  G_1(-,  g))  G 
Q(£l;g)  for  Q(Q;g)  =  {(q,z)  G  Y(Cl(g))  x  R(f 1(g))  :  b((q,  z),w;  g)  =  £(w;g),  Vw  G  V(Cl(g))}. 
We  in  addition  note  that  the  dual  energy  functional  admits  the  following  transformation:  for 
q  =  (det (T(g))T~T (g)q)  o  G(-;  g)  and  z  =  (det (T(g))z)  o  G(-;  g), 

Jd{(Q,z)‘,g)  =  -\  /  (q- D~1(g)q  +  c~1(g)zz)dx 

^  J 

q-  D~1(g)q  +  c~1zz^j  dx  =  Jd((q,z);  g), 

where  D(g)  and  c(g)  are  as  defined  for  a(-,  •;  g). 

We  summarize  the  computational  strategy.  We  first  recast  the  primal  and  dual  problems  on  Cl(g) 
as  those  on  the  reference  domain  Cl;  we  identify  the  associated  forms  a(-,  •;  g),  £(■;  g),  Jp{-;  g),  &(•,  •), 
and  Jd{-\g)-  We  then  apply  the  method  developed  in  Section  2  (with  the  extensions  introduced  in 
Sections  3.1  and  3.3)  in  the  reference  domain  to  construct  the  bounds  J^{g)  and  J^(g). 
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3. 5.  Extension  5:  Linear  Elasticity 

We  now  consider  an  extension  of  the  method  to  linear  elasticity.  For  clarity  in  this  section,  we 
employ  the  index  notation  with  the  implied  summation  on  repeated  indices.  We  first  introduce  a 
space  of  vector- valued  functions 

v  =  {v  G  ( H\n))d  :  Vi\rD,  =  0,  i  =  1, . . . ,  d}, 


where  T n.,  is  the  boundary  with  the  homogeneous  Dirichlet  condition  on  the  i-tli  equation.  The 
bilinear  form  is 

a(w,v;n)  =  /  eki(v)Ckuj(fj,)eij(w)dx-, 

Jq 

here  e(v)  G  ( L2(Ll))dxd  is  the  strain  tensor  given  by  €ij{v)  =  +  g^1),  and  C{ji)  is  the  rank-four 

stiffness  tensor.  The  linear  form  is 


t(w)  =  /  fkVkdx  +  /  gkvkdx, 


where  /  G  (L2(Q))d  specifies  the  body  force,  and  g  G  ®f=1L2(Pjv,i)  specifies  the  traction  on 
boundaries.  The  associated  energy  functional  is  Jp(w\  g)  =  ^a(w,w,g)  —  £(w).  Again,  owing  to 
the  variational  structure  of  the  problem,  the  construction  of  an  upper  bound  is  straightforward. 

To  construct  a  lower  bound,  we  first  introduce  a  broken  space  of  vector-valued  functions  V 
associated  with  V  and  a  space  of  rank- two  tensor- valued  functions  Y  =  ( V)d .  We  then  redefine  the 
bilinear  form 


b(q,w) 


/  qkitki{w)dx ,  Vq  G  Y,  \/w  G  V  . 

Jn 


We  next  redefine  the  dual  feasible  space 


Q  =  {q  G  Y  :  b(q,w )  =  £(w),  \/w  G  V}. 


We  finally  redefine  the  dual  energy  functional 


Jd(q ;  ^  =  ~\  fn  VkiCkUj  iv)qijdx,  VqeY  . 

We  may  then  show  that  J{n)  >  Jd{q\g)i  Vg  G  Q,  following  a  vectorized  version  of  the  proof  of 
Proposition  1. 

In  order  to  construct  a  finite  element  approximation  of  the  dual  variable,  we  follow  the  procedure 
of  Pares  et  al.  [10] .  We  first  extend  the  definition  of  the  trace  space  A  to  vector- valued  functions 
(in  the  same  manner  as  V).  We  accordingly  refine  the  bilinear  forms  c(-,  •)  and  &(•,  •)  as 


c(w,  x)  =  V  /  crKwkXkds, 
KdTh  JdK 

b(q,w)=  ^  /  qkiCki(w)dx , 


\/w  G  V,  Vx  G  A 
Vg  G  Y,  \/w  G  V  . 
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We  may  readily  compute  the  vector- valued  equilibrating  inter  elemental  flux,  G  A,  as  defined 

by  (7).  On  the  other  hand,  the  construction  of  the  dual  feasible  functions  via  the  localization  is 
more  complicated  than  for  the  diffusion  equation,  not  because  it  is  a  vector  equation,  but  because 
the  map  from  the  space  of  gradients  to  the  strain  is  not  bijective.  Nevertheless,  a  piecewise  Pp 
polynomial  representation  of  a  function  in  Q  over  each  may  be  found,  for  any  f\K  G  (Ppp  (n))d, 

g |7  G  (Pp/ (7))^,  and  A|7  G  (Pp(7))d,  Vy  G  dn,  using  a  subgrid  computation  procedure  described  by 
Pares  et  al.  [10];  we  omit  the  presentation  of  the  procedure  for  brevity.  (See  also  Pled  et  al.  [12] 
for  alternative  strategies.) 

While  the  finite  element  approximation  of  the  dual  variable  is  more  complicated,  the  reduced 
basis  approximation  requires  only  minor  modifications  from  the  diffusion  equation  case  considered 
in  Section  2.  This  is  because,  for  =  Yln=i  Yln'=i  (^(n'^^n'nPNn,  the  condition  for  pjy  G  Q  is 
En=i  Z£=1  Un’nPNn  =  1  as  before. 

In  the  offline  stage,  we  first  compute  snapshots  pft (pin1))  £  Q-  n>  =  1,...,-/V.  We  then  or- 
thonormalize  the  functions  to  form  the  basis  {vtf}n=\  for  the  reduced  basis  space  Yn  and  the 
change  of  basis  matrix  cu  G  WNxN .  We  next  sum  the  leading  index  of  the  change  of  basis  matrix 
to  form  w®um  =  Yhn’= \urin-  We  then  compute  parameter-independent  matrices  Kyu  G  M.NxN, 
i,j,k,l  =  1, . . . ,  d,  with  entries 

(K kilj)mn  =  /  Vm,ki1Tn,ljd,X  ’ 

Jn 

here  rjm^i  denotes  the  (k,  i )  tensor-entry  of  the  m-th  basis  function.  In  the  online  stage,  we  first 
form  a  parameter-dependent  matrix  K(/i)  G  M.NxN  defined  by 


K(/x)  =  Kkilj- 

k,i,l,j= 1 

we  then  solve  a  quadratic  program 

Jn(p)=  sup  Jd(qN]l-i)  =  sup  Jd(qN-,g) 


Qn&Qn 


(cosum)Tl3N=l 


with  a  single  linear  constraint  and  form  the  lower  bound  as  before.  The  computational  complexity 
in  the  online  stage  is  0(d4N 2)  +  0((N  +  l)3). 

3.6.  Extension  6:  Affine  Geometry  Transformation  (Linear  Elasticity) 

We  now  consider  an  application  of  the  method  to  linear  elasticity  problems  with  affine  geometry 
mapping.  The  strategy  is  similar  to  that  for  the  reaction-diffusion  equation  presented  in  Section  3.4; 
however,  we  must  now  transform  the  vector-valued  primal  variable  and  rank-two  tensor-valued  dual 
variable  in  a  consistent  manner.  As  before,  the  transformation  from  the  parameter-independent 
reference  domain  tt  to  the  parameter-dependent  transformed  domain  Q(/j)  is  denoted  by  x  = 
G(x;  n)  =  T(fi)x  +  xq(i u).  To  avoid  the  notational  clutter,  we  simply  state  T(/i)  as  T  from  hereon. 


We  first  note  the  transformation  of  the  strain  tensor 
w)  where  ekl(w)  =  \ 

wk  =  ( Tikwi )  o  G(-;  n)  and  vk  =  ( Tikvi )  o  Gf;  //), 


TrlTrl 


for  wk  =  (Tikwi)  o  €ij(w )  = 

+  j  •  We  next  transform  the  bilinear  form  for 


a(w,v;p)=  /  eki(v)Ckuj(n)eij(w)dx  =  =  a(w,  v\ p), 
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where  CW -,(//)  =  det (T ) Tgkl TrJ Ta  1  TrJ Csrntn{lJ) ■  We  then  transform  the  linear  form:  for  vk  = 
(Tikvi)  o  G(-\n), 


£{w) 


fkvkdx  + 


gkVkds 


fkVkdx  +  gkvkds  = 
J  £1  J  T  jv 


where  fk{p)  =  (det(T)Tfci1/i(/x))oG(-;  At)  and  gk(n)  =  (det(TrJV)Tw15i(/x))oG(-;  m).  We  in  addition 
define  the  primal  energy  functional  on  the  reference  domain:  Jvf\  /j)  =  ja(-,  •;  //)  —  £(•;  /j). 

We  now  consider  transformations  associated  with  the  lower  bound  construction.  We  transform 
the  bilinear  form  which  again  is  parameter  dependent  due  to  the  presence  of  as 

follows:  for  qki  =  (det(T)Tfct1T^1g^)  o  Gf;/j) 


b(q,  w,  n)  =  /  qki€ki{w)dx  =  qkieki(w)dx  =  b(q,w). 

jQ(fl) 

We  again  note  that  the  transformed  bilinear  form  &(•,  •)  is  independent  of  the  parameter  fi.  We 
then  introduce  the  space  of  dual  feasible  functions  on  the  reference  domain  Q(Cl\  n)  =  {q  G 
Y"(f2)  :  b(q,w)  =  £(w',n),  \/w  G  We  readily  confirm  that,  if  q(n)  G  Q(f2;/x),  then 

Qki  =  ((det(T))-1rfciTygy)  o  G_1(-,/r)  G  for  =  {?  £  Y(Q(n))  :  b(q,w\n)  = 

£(w,n),  \/w  G  V(f2(/i))}.  Finally,  we  note  that  the  dual  energy  functional  admits  the  following 
transformation:  for  =  (det(T)Tfc£1T/y1gi:))  o  G(-;//), 

=  [  QkiGk^Afi)qijdx  =  I [qkiCj^Aii)qijdx  =  Jd(q-,iJ ), 

where  C1  is  as  defined  for  5(-,  ■;  /x). 

The  computational  strategy  follows  that  of  the  reaction-diffusion  case  in  Section  3.4.  We  first 
recast  the  primal  and  dual  problems  on  54  (/i)  to  those  on  the  reference  domain  Q;  we  identify  the 
associated  forms  a(-,  •;  /r),  £(•;  //),  /r),  6(-,  •),  and  J7d(-;  g).  We  then  apply  the  method  developed 

in  Section  2  (with  the  extensions  introduced  in  Sections  3.1  and  3.5)  in  the  reference  domain  Q  to 
construct  the  bounds  ./^(/i)  and  J^(n). 


4.  Results 

4-1.  One- Dimensional  Reaction- Diffusion  Equation 

We  first  apply  the  exact  certification  technique  to  a  parametrized  reaction-diffusion  equation: 
—/lAu  +  u  =  x  over  fl  =  (0, 1)  with  homogeneous  Dirichlet  boundary  conditions  at  x  =  0  and 
x  =  1.  The  diffusion  coefficient  takes  on  fi  G  V  =  [10— 3 , 1] .  Our  finite  element  spaces  consist 
of  ne  equal-sized  P2  elements  (A f  =  2 ne  —  1).  Our  reduced  basis  space  is  constructed  from  the 
finite  element  solution  evaluated  at  N  Chebyshev-Lobatto  nodes  of  V  mapped  by  a  logarithmic 
transformation  [9].  The  exact  energy  J(g)  is  computed  using  a  P30  pseudo-spectral  discretization. 
Our  reduced  basis  formulation  exercises  the  extension  considered  in  Sections  3.3. 

The  bound  behavior  for  the  ne  =  128,  N  =  3  case  is  shown  in  Figure  1.  The  P2  finite  element 
discretization  with  ne  =  128  elements  is  sufficient  to  obtain  the  maximum  error  over  the  parameter 
range  of  O(10”8);  hence,  this  is  the  typically  assumed  reduced  basis  scenario  where  the  finite 
element  discretization  may  be  taken  as  the  “truth.”  As  a  result,  the  error  plot  over  the  parameter 


19 


(a)  energy  (b)  error 

Figure  1:  The  behavior  of  the  bounds  for  the  reaction-diffusion  problem  for  ne  =  128  and  N  =  3.  (J:  truth;  J^: 
reduced  basis  bounds;  J±,7Y  finite  element  bounds) 

domain  exhibits  the  typical  reduced  basis  method  error  behavior;  namely,  the  error  essentially 
vanishes  at  the  reduced  basis  parameter  evaluation  points  and  grows  away  from  the  points. 

The  bound  behavior  for  the  ne  =  4,  N  =  3  case  is  shown  in  Figure  2.  The  P2,  ne  =  4  finite 
element  discretization  provides  insufficient  resolution  particularly  for  a  small  /i;  this  is  reflected  in 
the  noticeable  finite  element  bound  gap.  The  quality  of  the  reduced  basis  bound  is  limited  by  the 
inadequacy  of  the  underlying  finite  element  discretization. 

Table  1  summarizes  the  maximum  error  of  the  reduced  basis  bound  over  S  C  T>,  where  the 
surrogate  subset  H  consists  of  1001  points  equidistributed  over  V  in  the  logarithmic  sense.  The 
error  in  the  bound  is  expressed  as  a  function  of  the  physical  space  resolution  (as  reflected  in  ne) 
and  the  parameter  space  resolution  (as  reflected  in  N).  We  first  report  that  the  upper  and  lower 
bounds  are  indeed  rigorous  bounds  of  the  exact  energy  J(/r):  J^)(/i)  —  J(/i)  >  0  and  J  —  >  0 

for  all  [i  E  T>,  ne,  and  N.  We  next  comment  on  the  two  extreme  cases  shown  in  the  table:  for 
ne  =  4,  the  physical  space  resolution  limits  the  bound  quality  independent  of  the  parameter  space 
resolution;  for  N  =  1,  the  parameter  space  resolution  limits  the  bound  quality  independent  of  the 
physical  space  resolution.  We  also  note  that  the  finite  element  bounds  converge  at  the  optimal 
rate  of  h2p  =  hA  for  this  smooth  problem.  Finally,  given  a  sufficient  finite  element  resolution  (for 
instance  for  ne  =  128),  the  reduced  basis  bounds  exhibit  an  exponential  convergence  with  N. 

4-2.  Planar  Linear  Elasticity 

We  now  consider  a  planar  (stress)  elasticity  problem  with  a  geometry  deformation;  we  refer 
to,  for  instance,  Rozza  et  al.  [15]  for  the  treatment  of  linear  elasticity  problems  in  the  standard 
reduced  basis  context.  We  consider  an  rectangular  elastic  beam  of  a  length  (normalized  with 
respect  to  the  height)  of  L/H  E  [2.0,  8.0]  composed  of  a  material  with  a  (non-dimensionalized) 
Young’s  modulus  1.0  and  Poisson’s  ratio  of  0.3.  The  beam  is  fully  clamped  on  one  end,  and  we 
apply  an  unit  tangential  traction  on  the  other  end;  all  other  boundaries  are  traction  free.  We 
then  compute  upper  and  lower  bounds  of  the  elastic  energy  as  the  length  is  varied;  recall  that, 
by  Proposition  3,  the  energy  bound  gap  may  be  used  to  construct  a  bound  of  the  energy  norm 
of  the  field.  Our  finite  element  spaces  consist  of  P3  elements  on  a  sequence  of  uniformly  refined 
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(a)  energy  (b)  error 


Figure  2:  The  behavior  of  the  bounds  for  the  reaction-diffusion  problem  for  ne  =  4  and  N  =  3.  (J:  truth; 
reduced  basis  bounds;  finite  element  bounds) 
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Table  1:  Convergence  of  the  reduced  basis  upper  and  lower  energy  bounds  (and  the  finite  element  bounds)  for  the 
reaction-diffusion  problem  with  the  number  of  elements  in  the  underlying  P2  finite  element  discretization  ( ne )  and 
the  reduced  basis  space  dimension  (N). 
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Table  2:  Convergence  of  the  (normalized)  reduced  basis  bound  gap  max^gac cd(T^(/t)  —  J^(/i))/|  Jref(/i)|  for  the 
linear  elasticity  problem  with  the  number  of  degrees  of  freedom  of  the  underlying  P3  finite  element  discretization 
(A f)  and  the  reduced  basis  space  dimension  (N).  The  reference  value  for  normalization  (and  not  assessment  per  se), 
|  Jref(p)\,  is  computed  using  a  Af  =  591,745  finite  element  solution. 


uniform  meshes.  Our  reduced  basis  space  is  constructed  from  the  finite  element  solution  evaluated 
at  N  Chebyshev-Lobatto  nodes  of  V  mapped  by  a  logarithmic  transformation.  The  reduced  basis 
formulation  exercises  the  extensions  developed  in  Sections  3.1,  3.5,  and  3.6.  The  number  of  terms 
in  the  affine  expansion  of  the  right-hand  side  (11)  is  S  =  1;  hence  we  may  readily  find  the  solution 
for  any  N  >  1. 

Unlike  the  reaction-diffusion  problem  considered  in  Section  4.1,  the  exact  solution  of  this  planar 
elasticity  problem  is  not  smooth;  in  particular,  in  the  presence  of  the  90°  corners  with  the  fully- 
clamped  and  traction-free  interfaces,  the  solution  can  be  shown  to  be  in  (i73/2+e(fl))2,  e  >  0,  for 
t  €  (L2(D))2  [14].  Hence,  we  will  not  attempt  to  compute  the  exact  output  for  this  case;  we  instead 
report  the  convergence  of  the  bound  gap,  —  ./() (ji) ,  with  the  finite  element  resolution  A f  and 

the  reduced  basis  dimension  N.  Note  that  this  is  in  fact  the  mode  of  operation  for  many  complicated 
practical  problems  and  is  precisely  when  the  rigorous  bounds  for  infinite-dimensional  variational 
problem  is  important:  the  computation  of  the  “exact”  solution  is  unreasonably  expensive  and  yet 
we  wish  to  have  guaranteed  certainty  in  our  computation. 

Table  2  shows  the  convergence  of  the  maximum  (normalized)  bound  gap  over  E  C  V,  where 
the  surrogate  subset  E  consists  of  201  points  equidistributed  over  V  in  the  logarithmic  sense.  The 
bound  gap  is  expressed  a  function  of  the  physical  space  resolution  (as  reflected  in  AT)  and  the 
parameter  space  resolution  (as  reflected  in  N ).  For  N  =  1,  the  parameter  space  resolution  limits 
the  bound  quality  independent  of  the  physical  space  resolution;  for  M  =  175,  the  bound  quality 
is  largely  limited  by  the  finite  element  resolution.  We  again  observe,  for  A f  =  148,417,  a  rapid 
initial  convergence  of  the  bound  gap  with  N.  Due  to  the  spatial  irregularity  mentioned  above,  the 
convergence  of  the  bound  with  the  uniform  mesh  refinement  is  rather  slow. 

5.  Summary  and  Discussions 

We  propose  a  reduced  basis  certification  strategy  that  provides  rigorous  upper  and  lower  bounds 
of  the  energy  associated  with  the  exact  infinite-dimensional  weak  solution  of  parametrized  steady 
symmetric  coercive  equations.  The  upper  bound  construction  is  based  on  the  usual  variational 
argument  over  a  subspace;  the  lower  bound  construction  is  based  on  a  reduced  basis  approximation 
of  the  dual  variable  that  satisfies  exactly  the  dual  feasibility  conditions.  We  then  consider  various 
extensions  of  the  basic  technology.  The  formulation  yields  truly  rigorous  certificates  of  the  reduced 
basis  energy  prediction  (and  the  energy  norm  of  the  solution)  without  any  assumptions  as  regard 
the  accuracy  of  the  underlying  finite  element  discretization. 
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We  identify  several  future  research  directions  that  address  some  of  the  limitations  noted  in  the 
Introduction.  We  also  reference  relevant  past  work  developed  in  different  contexts.  First  is  the 
extension  to  coercive  but  nonsynnnetric  equations  and  to  other  functional  output  quantities;  we 
refer  to  Sauer-Budge  and  Peraire  [17]  for  the  formulation  in  the  finite  element  context.  Second  is 
the  extension  to  noncoercive  equations;  we  refer  to  the  work  of  Peraire  and  Patera  [11]  for  a  related 
formulation  which  provides  asymptotic  bounds  in  the  finite  element  context.  Third  is  the  extension 
to  nonlinear  equations;  for  a  limited  class  of  nonlinearities,  it  may  be  possible  to  incorporate  the 
technique  developed  in  Machiels  et  al.  [8]  in  the  finite  element  context.  Fourth  is  the  extension  to 
nonaffine  parameter  dependence;  we  may  incorporate  the  empirical  interpolation  method  [2]  and 
the  associated  error  bounds  [3].  Fifth  is  the  simplification  of  the  finite  element  equilibration  pro¬ 
cedure;  we  may  incorporate  certain  discontinuous  Galerkin  methods  which  automatically  produce 
equilibrated  fluxes  [19].  Sixth  is  the  application  to  large-scale  engineering  systems;  we  may  develop 
an  approach  similar  to  the  component-based  reduced-basis  method  [4],  Although  all  these  exten¬ 
sions  should  in  principle  be  possible,  none  is  simple  and  furthermore  the  complementary  energy 
approach  is  clearly  most  advantageous  in  the  symmetric  coercive  context. 

We  also  identify  an  important  related  development  which  is  more  widely  applicable:  an  adap¬ 
tation  strategy  that  controls  both  the  physical  space  error  due  to  the  lack  of  the  finite-element 
resolution  and  parameter  space  error  due  to  the  lack  of  reduced-basis  resolution.  We  note  two 
recent  work  that  demonstrate  the  effectiveness  of  spatial  adaptivity  in  the  reduced  basis  context: 
the  work  of  Steih  and  Urban  [18]  which  uses  an  adaptive  wavelet  discretization;  our  work  which 
uses  a  minimum- residual  mixed  method  [20].  The  latter  formulation  is  built  on  duality  but  not  on 
complementary  energy.  In  this  adaptive  context,  we  can  consider  a  larger  class  of  equations,  as  the 
emphasis  is  on  the  development  of  an  optimal  approximation  even  if  we  cannot  confirm  accuracy 
with  complete  rigor. 
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